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ABSTRACT 

The quest for the nature of interatomic interactions in solids is of 
paramount importance as it leads to an understanding of their 
thermodynamical, elastic and numerous other physical 
properties. The thermodynamical properties either completely 
dictate the nature of the response of the materials or control the 
driving force for kinetic steps. The proposal presents an 
introduction to quasi harmonic Debye model and its applicability 
aspect in the understanding of thermodynamical, elastic and 
numerous other physical properties of wide range of materials. 



I. INTRODUCTION 

The equation of state (EOS) and the 
chemical potential are two of the key 
thermodynamic properties of a solid. The 
EOS of a given crystalline phase determines 
its behavior with respect to changes in the 
macroscopic variables, mainly pressure (p) 
and temperature (T). Experimental 
measurements are daily performed to 
determine the EOS of new materials or to 
extend it to wider p and/or T ranges. The 
results are usually expressed in terms of a 
reduced number of parameters by fitting the 
pressure-volume experimental data to an 
empirical equation. The chemical potential p. 
(equivalent to the molar Gibbs function, G m ) 
is the magnitude governing phase stability 
and change. Although its measurements are 
not so common, its effects in terms of phase 
transition properties and ranges of phase 
stability are some of the most important in 
state-of-the-art experimental solid state 



science. On the theoretical side, the 
determination of the EOS and the chemical 
potential from first principles are also two of 
the main objectives of the Physics and 
Chemistry of crystals. To obtain them, one 
has to pay particular attention to the concept 
of thermodynamic equilibrium state. 
According to standard thermodynamics, if 
the system is held at a fixed T and suffers a 
constant and hydrostatic p, the equilibrium 
state is the one that minimizes the 
availability or non-equilibrium Gibbs 
energy of the crystal phase 1 , 

G*(x; p, T) = E(x) + pV(x) + A vib (x; T) 
with respect to all internal configuration 
parameters. These configuration parameters, 
gathered in the configuration vector x, 
include all the relevant geometric 
information for the given crystal structure, 
i.e. the independent unit cell lengths and 
angles of this phase, and all the free 
crystallographic coordinates of the atoms in 
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non-fixed Wyckoff positions. On the right- 
hand side of above equation, E(x) is the total 
(or cohesive, if we set the energy zero at the 
infinitely separated components) energy of 
the crystal, which corresponds to the 
potential energy surface as determined by 
most electronic structure or atomistic 
calculations. The second term, pV, 
corresponds to the constant hydrostatic 
pressure condition. Finally, the third term, 
A vib , is the vibrational Helmholtz free 
energy, which includes both the vibrational 
contribution to the internal energy and the 
-T S constant temperature condition term 
(since a perfect crystal is considered whose 
only internal degrees of freedom are 
vibrations, S=S v ib)- The rigorous statistical 
calculation of A vib requires knowledge of the 
exact vibrational levels, but it is customary 
to introduce the quasi-harmonic 
approximation 1 , 

A V ib(x; T) = 

[-hte + kTln(l - ekT)]g( x; o))dco 
J 0 2 

where g(x;ro) is the phonon or vibrational 
density of states. The term quasi-harmonic, 
as different from the rigid harmonic 
approximation, remarks that the density of 
states is allowed to vary with the crystal 
configuration, thus including anharmonic 
contributions to a certain extent. 

Summarizing the above paragraph, a 
calculation of (p,T) equilibrium states 
requires knowledge of E(x), V (x), and 
g(x;co). Since V (x) can be readily computed 
for any given phase, and the potential energy 
surface is the main product of solid state 
computational codes, a first approach in 
which the vibrational effects (both finite- 
temperature and zero-point) are neglected is 
commonly used. In this static approach, the 



static non-equilibrium Gibbs energy, 
G*static(x;p) = E(x) + pV(x), is minimized 
with respect to x. Although the full 
minimization can usually be attained, in 
cases of high dimensionality and when using 
ab initio methods a restricted minimization 
can also be employed, for example fixing 
some of the internal parameters, cell angles, 
or ratios of cell lengths. Once G*static has 
been fully or partially minimized, the static 
EOS of the solid is recovered as V (p) = V 
(x opt (p)), and its chemical potential as G(p) = 
G* st atic(x 0 pt(p);p)- Although these results are 
quite valuable and give good insights into 
the general behavior of crystalline phases, 
they are not comparable to the experimental 
results. Even for low T, thermal effects are 
present both in the form of zero point and 
finite temperature contributions, and their 
inclusion in the theoretical study of crystals 
is very relevant, not only when studying 
temperature dependences but also for fixed- 
T studies. Some methods have recently been 
developed to fully minimize G*(x;p,T) with 
respect to x within the quasi harmonic 
approximation when simple pair potentials 
are used in the crystal simulation 2(a) . 
However, when a first principles method is 
used, the full minimization of G*(x;p,T) is 
currently infeasible. This fact has motivated 
M.A. Blanco et a/. 2(b) to develop a non- 
empirical model to include thermal effects in 
a general way after a possibly expensive 
calculation of the potential energy surface of 
a crystalline phase 3 " 7 . M. A. Blanco et al. 2(b) 
have thus made it work with a minimal 
amount of data, in this case a set of 
computed pairs of energies and volumes per 
unit formula, that correspond to the lowest 
energy configuration x opt compatible with 
that volume. The simplicity of the 
requirements makes it completely general, 
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although it is expected to work better for 
highly symmetric crystals. 

II. THE QUASI-HARMONIC DEBYE 
MODEL 

Basic input of the model: static 
optimizations 

Let us assume that an appropriate 
(ab initio 812 or otherwise) approximate 
algorithm to compute E(x) is available. 
Then, before the quasi-harmonic Debye 
model can be properly applied, this 
multivariable surface has to be transformed 
into an E(V) curve. This can be done in two 
equivalent ways. On the one hand, since the 
volume is usually a quite simple function of 
the lattice parameters, one can minimize 
E(x) for a set of fixed volumes, for example, 
optimizing the internal coordinates, angles, 
and cell length ratios (c/a and b/a); writing V 
= abcf(a,p,y)/Z = a 3 g(a, P, y, b/a, c/a), where 
Z is the number of molecules per unit cell, 
and f and g depend on the crystal system and 
Bravais lattice of the phase, the last cell 
length is then univocally defined for each 
volume to be a = [V/g(a, p, y, b/a, c/a)] 1/3 . 
This produces a curve E(V) = E(x opt (V)) that 
contains the points of the surface which have 
lowest energy for different constant V 
hyperplanes. On the other hand, one can 
perform static optimizations at several 
constant pressures, using the G* static (x;p) = 
H* stat i c (x;p) = E(x) + pV(x) function to 
obtain x opt (p). As shown above in the 
introduction, this leads to the static EOS V 
(p) = V (x opt (p)), and thus M. A. Blanco et 
a/. 2(b) obtained E(p) = E(x opt (p)). Since the 
EOS was a table of (V,p) values, they invert 
it to obtain (E,V) pairs, as they wanted. Both 
procedures are equivalent, since H* static is 
nothing but the Lagrange function for the 



minimization of E at constant V, being p the 
corresponding Lagrange multiplier. Thus, 
both E(V) curves are one and the same, and 
choosing one or the other of these 
procedures will only change the actual 
points selected from the single curve. Since, 
in what follows, the E(V) curve is the key 
input data for several numerical 
minimization procedures, V points with a 
homogeneous distribution should be used, 
and they should span values both above and 
below the static zero-pressure equilibrium 
volume V 0 . Thus, if no previous information 
on the V (p) EOS is available, fixed volume 
optimizations are preferred, selecting an 
appropriate volume grid. This avoids the 
selection of a pressure grid that should 
include unphysical negative pressures, in 
order to sample points above V 0 . An 
important word of caution should be given 
here. Since by selecting an univocal 
correspondence x opt (V) between configure - 
ations and volumes, given by the static 
calculation, M.A. Blanco et al. 2(h) effectively 
assumed that thermal effects over the 
configuration can be reduced to the thermal 
effects over volume. Having x = x opt (V ) is 
the case in the static calculation, and usually 
holds for pressure -induced effects at 
constant temperature, too; it can even be a 
good approximation for temperature -induced 
effects, but this will only be true when the 
vibrational effects act as a hydrostatic term, 
i.e. are isotropic. 

Thermal effects: quasi-harmonic Debye 
model 

After the static calculation, the non 
equilibrium Gibbs function G*(V;P,T) can 
be written in the form of 5 ' 2(b) 
G*(V; P, T) = E(V) + PV + A vib (0(V); T), 
where E(V) is the total energy per unit cell, 
PV corresponds to the constant hydrostatic 
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pressure condition, and A V ib (0(V);T) is the 
vibrational term, which can be written as 

A vi b(©(V);T) = 

nKT X [-® + 3 In (l - e^) - D (0 / T) ] 

where n is the number of atoms in the 
molecule, and the Debye integral D(0/T) is 
defined as 2(b) 



a = 



3 ry x 3 

y 3 J 0 e x - 1 



dx. 



For an isotropic solid, 0 is expressed by 

0 = r[67r 2 V 1 /2n] 1 / 3 f( o ) 

where M is the molecular mass per formula 
unit, Bs the adiabatic bulk modulus, which 
can be approximated by the static 
compressibility 2(b) 

d 2 E(V) 

B s~B(V) = V(-^) 



and the Poisson ratio o and f(o) 5 ' ] are given 
by the following forms: 
3B-2G 

o = 



6B + 2G' 
[f(o)] 3 = 3[2(- 



2 l + o 3/ 



3 1 - 2a J 

11 + 03/ 1 

Therefore, the non-equilibrium Gibbs 
function G*(V;P,T) as a function of (V;P,T) 
can be minimized with respect to volume: 

f^],T = 0 (1) 

The isothermal bulk modulus B T , the heat 

capacity C v and the thermal expansion 

coefficient (a) are expressed as 2(b) 

d 2 G*(V; P,T) 
B T (P,T) = V[ ^^]p,t 

C v = 3nk [4D (0/T) — ] 



yc v 

B T V 



where y is the Gruneisen parameter defined 
as 

Vln0(V) 
J ~~ dlnV 

III. COMPUTATIONAL SCHEME 

To simplify the process of 

minimization/derivation involved in the 

above equations, it is convenient to fit an 

appropriate analytical function of V to the 

numerical values of G*(V;T,p) and E(V). 

Since typical energy curves resemble 

anharmonically distorted parabolae, simple 

polynomials in V are commonly used. 

However, this possibility leads to large 

errors in the values of the fitted function, 

specially at high volumes. To minimize 

these errors, instead of V, M.A. Blanco 

et a/. 2(b) used the reduced length unit x = (V 

/V re f) 1/3 , where V re f is chosen as the value of 

V in the (E,V) table for which E is a 

minimum. Even with this improvement, the 

derivatives obtained were highly dependent 

on the polynomial order and on the end point 

values of the data set. To avoid this problem, 

they fit polynomials with a wide range of 

orders, and successively redo the fitting 

eliminating the endpoints from the input 

data. In this way, get a large number of 

different polynomials to estimate G*(V;T,p). 

Then, they averaged the fittings, assigning to 

the i th polynomial a weight equal to 

_ 2 _ 2 
pi = e "'/Die ffl i, where 

co i = [(ni + lVN^Sj/S,™] n ; is the 

polynomial degree, N ; the number of (E,x) 

data, 5i is the root mean square deviation in 

the i* fitting, and 8 min = minjSi}. In this 

way, they 2(b) obtained numerically stable 

values for both the minimum of the function 
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and its successive derivatives, while 
retaining the simplicity of the polynomial 
fittings. 

The above fitting strategy was 
applied both to the static (E,x) data and to 
the thermal (A*,x) data (A*(V;T) 
=E(V)+A vib (@(V);T)). It is worth 
mentioning that only the zero-pressure data 
at each T need to be fitted. Indeed, if f (x) 
=2j aiXi is the polynomial that best fits the 
(E,x) data, the best polynomial fitting for the 
(E + pV,x) data with p ^ 0 coincides with it 
after changing the third-order coefficient to 
a 3 "= a 3 +pV ref . This is evident when 
considering that pV = px 3 V ref . Once the E(x) 
function has an analytical form, 0(x) can be 
obtained by simple analytical derivation of 
the polynomial, and the (A*(x;T),x) data 
obtained. The fitting strategy then provided 
them 2(b) with the G*(x;p,T) polynomial 
function for any selected temperature and 
pressure. Then, solving Eq. (1) by a 
combination of bisection in the (G*, x) data 
and the Newton-Raphson method, since 
both the first derivative (the function to be 
zeroed) and the second derivative (the 
derivative of the function to be zeroed) are 
polynomials easy to obtain from the original 
G*(x) one. A further auxiliary polynomial fit 
was performed with the In© versus InV 
values, both as a means of accessing & at 
values other than those in the input, and as a 
means of obtaining the Griineisen parameter. 

IV. EMPIRICAL ANALYTICAL EOS 

The Gibbs program can fit analytical 
expressions to the p(V) data obtained by 
solving Eq. (1). These empirical EOS 
involve the zero-pressure bulk modulus, B 0 , 
and its first, second, . . . pressure derivatives 
(B 0 \ Bo"", . . .) as fitting parameters, thus 
providing an alternative derivation of the 



compressibility. These analytical express- 
ions can easily be manipulated to derive 
most thermodynamic properties at any T and 
p. The static analytical EOS can be likewise 
used to compute 0(V) after a corresponding 
fit to the static p(V) values. As a means of 
estimating the quality of the EOS fitting, the 
explicit p(V) expressions can be analytically 
integrated to obtain the corresponding 
energy (E in the static calculations, or the 
Helmholtz non-equilibrium energy function 
A* in the thermal calculations), which can 
then be directly compared with the original 
energy volume data. Presently, three 
possibilities are incorporated in Gibbs: the 
Vinet et al. EOS 14 , the Birch-Murnaghan 
EOS 13 , and the spinodal EOS 15 . 

The Vinet EOS 

This EOS connects the p(V) data through the 
relation 

ln[px 2 /3(l - x)] = In B 0 + a(l - x), 

x= (V/Vo) 1/3 

where V Q = V(0,T) is the zero pressure 
equilibrium volume and In Bo and a=3(Bo"- 
l)/2 are the fitting parameters. 
B T = -x- 2 B 0 e a(1 - x) f(x) 

Bf = (^-) T = l/3[(ax+2) 

-xf (x)/f(x)], and B T " = (^) T = 
x/9B T{ x^(x)/ f(x) _ a + r(x)/ f(x)[1 _ 

x ^ w here f(x) =x-2-ax(l-x). 

The Birch Murnaghan EOS 

In the second order Birch 
Murnaghan EOS the p(V) data are connected 
by the equation 
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F(f) = 



3f(l+2f) 5 / 2 



- V2 



Zi=o a i f' 



where f=l/2(x -1) (x having the same 
meaning as in the Vinet EOS), and a Q , a 1; 
and a 2 are fitting parameters. The 
expressions for B T , B-f and B T ~ at zero 
pressure values are given by 
B 0 = a 0 

B x 0 = 4 + (2al/3B 0 ) 

r2a 2 /143\1 



/B 0 



The spinodal EOS 



The spinodal EOS proposed by Baonza et 
al. 15 connects the (p,V) data through the 
relation 

V L = V sp exp {-[JL^](p-p sp )i-P} 
where 

Vs p =v ° exp{ (T^w; 

In the above equation, p sp and v sp are the 
spinodal pressure and volume respectively. 
The fitting parameters are p sp , K*, and p (it 

is customary to use P=0.85). In terms of 
them, we have 

B 0 = [K1- 1 (-Ps P ) P 
B 0 ^ = (-Psp) -1 PB 0 . 

At any T and p, B T , B T " and B T "" are 
straightforwardly given by 

B T =[K*]- 1 (p-p sp ) P 
Bf = P[K*]- 1 (p-Psp) P " 1 

B T " = p(p-i)[r]- 1 (p-p sp ) p - 2 

The expression of B T is the one used to 
optimize (-p sp ), K*, and p from the B T 
values. Putting it in the form 



lnB T = pin(p-p sp )-lnK* 

And given a trial value of (-p sp ) and the set 
of data (p,B T ), K*and p are obtained solving 
the corresponding linear least squares 
problem. The optimum (-p sp ) is trivially 
found by a monodimensional minimization. 

V. OTHER THERMODYNAMIC 
PROPERTIES 

After the equilibrium state for a 
given (p,T) has been obtained, other 
thermodynamic properties can also be 
evaluated by using the corresponding 
equilibrium volume in the appropriate 
thermodynamic expressions. For instance, 
the vibrational internal energy (U v ;b), heat 
capacity (C v>vib ), and entropy (S vib ) in the 
quasi-harmonic Debye model are given 
by 2 *) 

9 0 r@/ 
Uvib = nKT[--+3D^ U / T) ] 



C v , vib = 3nk [4D (0 / T , - \r^- 
I J e /t _ i 

S vib = nk [4D (0 / T) - 3 In (l - e" 8 ^)] 

Another relevant property is the Gruneisen 
parameter, defined as 



Y = -■ 



Vln0(V) 
dlnV 



The only explicit dependence of this 
parameter is on V. However, since the 
derivative must be evaluated at the 
equilibrium volume at each T and p, it has 
an implicit dependence on these two 
variables. Although the above equation can 
be used directly to obtain y, it is more 
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accurate to derive it from the Mie-Gruneisen 
equation 

U vib 

P-Pstatic = y~^~ 

The second member of this equation 
represents the thermal contribution to the 
pressure. Finally, the thermal expansion (a), 
the heat capacity at constant pressure (C p>vib ), 
and the adiabatic bulk moduli are given 
by 2(b) : 

_ Y C v,vib 

a_ B T V 

Cp.vib = C v,vib(l + ayT) 
B s = B T (1 + avT) 
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